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We describe a numerical method for calculating the (3+1) 
dimensional general relativistic hydrodynamics of a coalescing 
neutron-star binary system. The relativistic field equations 
are solved at each time slice with a spatial 3-metric chosen to 
be conformally flat. Against this solution to the general rela- 
tivistic field equations the hydrodynamic variables and grav- 
itational radiation are allowed to respond. The gravitational 
radiation signal is derived via a multipole expansion of the 
metric perturbation to the hexadecapole {I = 4) order includ- 
ing both mass and current moments and a correction for the 
slow motion approximation. Using this expansion, the effect 
of gravitational radiation on the system evolution can also be 
recovered by introducing an acceleration term in the matter 
evolution. In the present work we illustrate the method by 
applying this model to evaluate various orbits of two neutron 
stars with a gravitational mass of 1.45 M© near the time of 
the final merger. We discuss the evidence that, for a realistic 
neutron star equation of state, general relativistic effects may 
cause the stars to individually collapse into black holes prior 
to merging. Also, the strong fields cause the last stable orbit 
to occur at a larger separation distance and lower frequency 
than previously estimated. 

PACS Numbers: 04.20.Jb, 04. 30. -fx, 47.75-ff„ 95.30.Lz, 
95.30.Sf 97.60.Jd, 97.80.-d 



I. INTRODUCTION 

Coalescing neutron stars are currently of interest for 
a number of reasons. Several neutron-star binaries are 
known to exist in the Galaxy (e.g. PSR 1913-f 16 H PSR 
2303-f46 §, PSR 2127+llC |], PSR 1534-f 11 g) whose 
orbits are observed to decay on a time scale of 1 — 3 x 10* 
yr. It has been recognized for some time ||5|-[l0| that 
the final stages of coalescence of such systems could be 
copious producers of gravitational radiation. This pos- 
sibility has recently received renewed interest with the 
development of next generation gravity-wave detectors 
such as cryogenic bar detectors |l^, the Caltech-MIT 
LIGO detector jl2j and its European counterparts, GEO 
and VIRGO (e.g. [Q) for which an event rate due to 
binary neutron star coalescence out to 200 Mpc could 
be >3 



per yr 16 1. It has also been proposed that 

such events (when integrated over the number of galax- 
ies out to high redshifts) could account for the observed 
event rate and energy requirements of gamma-ray bursts 
-|l9i. Coalescing neutron stars might even be signifi- 



cant contributors to heavy element nucleosynthesis in the 
Galaxy |^,|l|. 

For much of the evolution of a neutron star binary, the 
system should be amenable to a point source descrip- 
tion using post-Newtonian techniques |22-2j]. However, 
as the stars approach one another the gravitational fields 
become quite strong and hydrodynamic effects should be- 
come significant. Indeed, it is expected that the wave 
forms could become quite complex as the stars merge. 
This complexity, however, may be sensitive to various 
physical properties of the coalescing system [|l^ such as 
the neutron star equation of state. Hence, careful mod- 
eling is needed which includes both the nonlinear general 
relativistic effects and a realistic neutron star equation of 
state. Such calculations can be used as a foundation for 
extraction of the information contained in the detected 
gravity waves and as a framework in which to analyze 
possible gamma-ray burst models.. 

A computation of the hydrodynamic evolution is 
complicated, however, due to the inherently three- 
dimensional character of the orbiting system. To this end 
several attempts have been made to model the hydrody- 
namics of coalescence in either a L agrang ian smoothed- 
particle Newtonian approximation p5 26 1 or using con- 
ventional finite-difference methods in the post-Newtonian 
approximation ||2^-^^. It important to appreciate, how- 
ever, that as the two neutron stars coalesce the system 
becomes strongly relativistic, and the validity of New- 
tonian or post-Newtonian hydrodynamics may be ques- 
tionable. In the present work, therefore, we improve upon 
such calculations in that we explicitly include most of the 
effects of a fully general relativistic treatment. 

Some preliminary discussion of this work has been re- 
ported previously |31-3^. In this paper we present de- 
tailed discussion of our method of solving the relativis- 
tic field equations and hydrodynamics. As an illustra- 
tion and for comparison with existing calculations in the 
literature, we present orbit calculations for two neutron 
stars with a gravitational mass of 1.45 Mq each. We find 
that the last stable orbit occurs for a separation distance 
w 1.4 times larger and a frequency smaller than that esti- 
mated using the post-Newtonian approximation. We also 
find the surprising result that with a realistic equation 
of state, the strong fields may induce otherwise stable 
neutron stars to collapse into black holes prior to orbit 
instability and merging. 
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II. THE MODEL 



A. Coordinate System 



We start with the usual shcing of spacetime into a one- 
parameter family of hypersurfaces separated by differen- 
tial displacements in time-like coordinates as defined in 
the ADM or {3+1) formahsm 

For this work we have considered a number of possible 
3-space coordinate systems, e.g. polar, bipolar, spheri- 
cal, cylindrical. Ultimately, we have selected Cartesian 
X, y, z isotropic coordinates. This is a natural coordinate 
system for three dimensional problems, in that no special 
point or singularity is introduced. It thus avoids prob- 
lems associated with finite differencing near coordinate 
singularities. It also has the advantage that the relativis- 
tic field equations assume a simpler and more symmetric 
form. 

With this choice for coordinates, proper distance is 
expressed 



\c? - Pi(3')dt^ + 2l3,dx'dt + jijdx'dx^ (1) 



where the lapse function a is a multiplier which describes 
the differential lapse of proper time between two hy- 
persurfaces. In the Newtonian limit this quantity ap- 
proaches unity and is related to the Newtonian gravi- 
tational potential. The quantity (3i is the shift vector 
denoting the shift in space-like coordinates between hy- 
persurfaces. The quantity jjk is the metric of the 3- 
geometry. It specifies the distance between points within 
a hypersurface. 

Here we introduce an approximation that the 3- 
geometry is both conformal and flat. That is, we write, 



and 



7^J 



(2) 



(3) 



where the conformal factor is a positive scalar function 
describing the ratio between the scale of distance in the 
curved space relative to our flat space manifold, and 5jk is 
the Kronecker delta. This is an approximate gauge condi- 
tion which we will henceforth refer to as the conformally 
flat condition or CFC. This approximation is motivated 
both by the general observation that gravitational radi- 
ation in most systems studied so far is small 39 , and 
the fact that conformal flatness on each space-like slice 
considerably simplifies the solution to the field equations. 

To see the way in which the CFC allows us to solve the 
relativistic field equations, consider the exact equation 



(4) 



where Di is the three-space covariant derivative [ p7[ , and 
Kab is the extrinsic curvature describing the deformation 



of a figure as it is carried forward by one unit in proper 
time in a direction normal to a hypersurface. 

Equation fl) is well approximated by a conformal rep- 
resentation (g) only if the trace free part of the right hand 
side vanishes. Thus, a spatially flat 3-metric requires. 



(5) 



where we have also employed the maximal slicing condi- 
tion, tr{Kab) = as a gauge choice. 

We use Eq. (^) to determine the extrinsic curvature. 
A convenient consequence of this is that any geometry 
which is initially conformally flat, will remain confor- 
mally flat to the extent that energy in gravitational ra- 
diation is unimportant. Equation allows us to derive 
constraint equations for the lapse function and conformal 
factor as described in the next section. 

As a final condition, we take the coordinate system to 
be rotating in such a way as to minimize the matter mo- 
tion in the coordinate grid. This condition enhances the 
stability of the computation of the hydrodynamic evolu- 
tion. However, this is a nontrivial condition to impose 
in curved spacetime which we achieve by boundary con- 
ditions on /3* as described below. All relevant forces are 
computed first in nonrotating coordinates which are then 
transformed to update the the matter fields in a rotating 
grid. 



B. General Relativistic Field Equations 

For most gravitating systems stud- 
ied so far (e.g. |^,^), only a relatively small amount 
of energy is emitted by gravitational waves. Even for the 
merger of two black holes it is expected |^ that only a 
few tenths of a percent of the rest mass will be radiated 
away in gravitation. For the case of two neutron stars we 
would not expect any more radiation to be emitted dur- 
ing the last few orbits than for a two black-hole merger, 
i.e. during the inspiral, the radiated energy per orbit is a 
minuscule fraction of the energy in orbital motion. Fur- 
thermore, an explicit treatment of the radiation reaction 
is exceedingly difficult js^. To treat the effec ts of g ravi- 
tational waves we use a multipole formalism 
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44|. We 



use a radiation reaction potential in the hydrodynamics 
equations to account for the effect of gravity waves on 
the system. 

The implementation of this approximation means that, 
given a distribution of mass and momentum on some 
manifold, we first solve the constraint equations of gen- 
eral relativity (GR) at each time in the calculation for 
a fixed distribution of matter. Then we let the mat- 
ter and gravitational radiation respond to this geometry. 
That is, we evolve the hydrodynamic equations to the 
next time step under an assumption of "instantaneous 
gravity." However, at each time step we obtain a time 
symmetric solution to the field equations. 
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As an alternative to the explicit coupling of emitted 
gravitational radiation to the hydrodynamic and geomet- 
ric evolution of the system, the initial evolution of the 
system (while the gravitational radiation is a small per- 
turbation) can be approximated by stable orbits in the 
absence of energy and momentum loss due to gravita- 
tional radiation. One can then after the fact compute 
the expected gradual loss rate of energy and momentum 
in gravity waves. This latter approach is applied in some 
illustrative calculation presented here. In a future paper 
we will subsequently apply the former method to describe 
the late time merging and coalescence. 

Our basic approach to a solution to the GR equations 
will be to reduce all of the constraint equations to effec- 
tive flat-space elliptic equations which are amenable to 
standard techniques. In what follows we make use the 
usual natural units in which G = c = 1. Thus, at each 
time slice we can obtain a numerically valid static so- 
lution to the exact GR field equations and information 
on the hydrodynamic evolution and generation of gravi- 
tational radiation. However, the advance from one time 
slice to the next incorporates the approximation that the 
effects of gravitational radiation can be neglected. 



1. Hamiltonian Constraint 



We begin with the Hamiltonian constraint equation 
[37|. We use the forms of equations as given by Evans 
[38|. We show here that the Hamiltonian constraint and 
the maximal slicing condition {tr{K) = 0) can be com- 
bined so as to form elliptic equations for both the con- 
formal factor (j) and the product (at/)). 

The Hamiltonian constraint equation can be written. 



R — WirpH 



(6) 



where R is the Ricci scalar curvature, and 

PH^phW^-P , (7) 

where, p is the proper baryonic matter density, W is 
the generalization of the special relativistic gamma fac- 
tor {W = all* , where U^^ is the four- velocity) , P is the 
pressure, and h is the specific relativistic enthalpy. 



h = l + e + P/p , 



(8) 



with e the associated matter energy above the baryon 
rest energy, and P the matter pressure. 

The conformal scaling of the 3- metric, Eq. (^), defines 
a conformal metric and manifold (7, M) related to the 
physical metric and manifold (7,M) [see refs. |]37| , ^ ]. 
Covariant derivatives Di and Di on M and M can be re- 
lated by calculating the transformation of the Christoffel 
connections, 



- jk 



20- 



SID, 



l3kl 



(9) 



With this, the transformation of the Ricci scalar cur- 
vature is 



where R 



R = (i)-^R 
R{j) and R = 



(10) 



i?(7) and A = f'D.Dj. As 



mentioned in Section ( II A ), we choose a conformally flat 
metric, 7*-'' = S^^ , for which, f jj, ^ 0, £», ^ V, i? -> 

and A V^, the flat space Laplacian. 

Solving Eq. ( |l0| ) for A(/), and combining with the 
Hamiltonian constraint gives the desired form for an el- 
liptic equation for 0, 



WnpH 



(11) 



In order to put this constraint equation into a form 
which is useful for solution along with the hydrodynamic 
variables, we must introduce conformal scalings for the 
source terms. To do this, the equation of state is intro- 
duced through the adiabatic index F: 

P={V-l)pe , (12) 

where F is a function of state variables for each zone. 
With this equation of state, (0) becomes, 

(F-1)- 



pu = pW^ + peW 



TW 



W 



(13) 



For the hydrodynamic Lorentz contracted density, D = 
pW and energy E = peW, we introduce the following 
conformal scalings. 



E = (h-^^E 



(14) 
(15) 



The reasons for these choices will be clear when we con- 



sider the hydrodynamic equations given in section (H C) 
The extrinsic curvature is scaled by, 



which gives 



(16) 



(17) 



With the introduction of these scalings, the Hamilto- 
nian constraint can be written into the desired form. 



-2ttE{ TW 



1 



1 



W J" 8 

This can be written in a more familiar Poisson form, 

V20=-47rpi, (19) 

in which the source term can be identified in terms 
of physical hydrodynamic variables by transforming the 
conformal scalings in Eq, 



(18) 



DW + E[ tw- 



(r-1) 
w 



IGtt ^ 



(20) 
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2. Lapse Function 

We also use the Hamiltonian constraint together with 
the maximal slicing condition tr(k) = to obtain an 
elliptic constraint equation for the lapse function a. We 
begin with the identities, 

Aa = ^[<t>-\a(t>)] = D^D'[4>-\a(j))] (21) 

(22) 

Now in our conformally-flat metric one can write for any 
scalar function, and in particular for the quantity (acf)) 

A(a0) = 0-^A(a0) + 2(t>-'f^{D,^)iD,{a(f>)) . (23) 

Substituting this into equation ( |2^ ) gives, 

Aa = 0-^A(a0) + a(j}A{<j)-^) . (24) 

Now from the transformation properties of the Ricci cur- 
vature scalar (|To|), equation (|2j) can be rearranged as 



A{a(j}) = (t)^Aa + -acj)^ 
8 



R(t>- 



R 



(25) 



Rewriting the Hamiltonian constraint (|6|) to include 
the CFC and maximal slicing conditions, then leads to a 
flat-space elliptic equation in (a0). 



1 



WirpiSW^ - 2) 



16TTpe[3T{W^ + 1) - 5] + 7K,jK'^ 



In Poisson-like form this is 

V^{a(t)) = 47rp2 



(26) 



(27) 



with the source term written in terms of hydrodynamic 
variables as 



P2 



D{3W^ - 2) + E[3T{W^ + 1) 



W 



167r 



(28) 



A solution of equation ( p7| ) determines the lapse function 
after equation (|l^) is used to determine the conformal 
factor. 



The momentum constraints have the form [ |38| , 

D,{K'^ ~ f^K) ^ SttS^ . (29) 

Where Dj is the three-space covariant derivative ( pTf). 
5* is the contravariant material momentum density which 
is derived from the solution to the hydrodynamic equa- 
tions, Section (ilC). In our maximal-slicing conformally- 
flat conditions, the second term on the left hand side of 
Eq. ( p9| ) vanishes and we have, 

D,K'^ = 8ttS^ . (30) 

Using ( p^ and ^ it can be verified that 

D,K'' === (f)-^° b,{(j)^° K'i) . (31) 

Now converting our "conformally fiat condition" 
[i.e. Eq. (H)] from covariant derivatives to or ordinary 
derivatives, and inserting into equation (^ij) gives. 



-10 



-D, 



a 3 



Combining this with (pO[) then gives, 



(32) 



(33) 



where the source term is defined. 



= 167ra(/.'*S'^ + ^ 



3 



where ^ = a . 

Equation ( |33|) can be reduced to. 

Thus, by introducing a decomposition of /?' into 
the following two elliptic equations result 



(34) 



(35) 



(36) 



(37) 



3. Momentum Constraint 

With the lapse function and conformal factor deter- 
mined from the Hamiltonian constraint and maximal slic- 
ing condition, we then use the momentum constraints to 
find the shift vector. 



(38) 



These we use to determine the components of the shift 
vector. 

This is the desired result except for the fact that most 
of the momentum encompassed in equations (|34|) and 
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( P8| ) is simply the orbital motion of the binary. We there- 
fore wish to define a rotating coordinate system with a 
rotation-subtracted shift vector in which the non-orbital 
aspects of matter evolution and relativity (e. g. frame 
drag) can be more easily studied. To do this we decom- 
pose -B' into a frame-drag term, G* and an orbital motion 
term. 



(39) 



and subtract the orbital velocity of the coordinate system 
from both sides of equation (|35|). Ultimately, we write 
equation (BSh in the Poisson-like form, 



where 



P3 



Aa(l3^Si-0W{D + rE)] 



(40) 



(41) 



1 dS, d(3' 



d(3^ 
dx'- 



Since V^((jJ x r) = 0, B^ and G* can be used inter- 
changably in equation (|40|). 

As the meaning of orbital angular velocity becomes 
obscured in curved spacetime, uj in equation ( p9| ) takes 
on the meaning of a Lagrange multiplier which minimizes 
the matter motion with respect to the coordinate system. 
It only reduces to the orbital angular velocity in the New- 
tonian (i.e. r — !■ oo) limit. Confining orbital motion to 
the X, y plane, we determine the coordinate rotation fre- 
quency UJ at each time step from the weighted average of 
the matter four velocity and the frame-drag shift vector, 



/ dV{D + TE) 



a(xUy-yU^) 



f-<j>\xGy-yG^) 



JdV{D + TE){x^+y^) 



(42) 



This rotation is then subtracted from the velocities and 
added to the coordinate rotation, thereby maintaining a 
centering of the stars along y ^ 0. 

The fact that the constraint conditions on (111), (27), 



and (40) can be written in the form of flat-spaced Poisson 
equations, allows for these variables to be solved by fast 
numerical techniques as discussed below. However, their 
solution requires that boundary values for these variables 
be specified at distances relatively close to the neutron 
stars. Our method of deter mining the boundary values 
is described in section [I F 2 . 



C. Relativistic Hydrodynamics 

To solve for the fluid motions of the system in curved 
spacetime it is convenient to use an Eulcrian fluid descrip- 
tion We begin with the perfect fluid stress-energy 
tensor, which in covariant form can be written. 



T^. = ip + pe + P)U,,U, + Pg^, 



(43) 



By introducing a set of Lorentz contracted state vari- 
ables it is possible to write the relativistic hydrodynamic 
equations in a form which is reminiscent of their New- 
tonian counterparts. The hydrodynamic state variables 
are: The coordinate covariant baryon mass density, 



D = Wp , 
internal energy density, 

E = Wpe , 

three velocity. 



(3^ 



and the momentum density, 

S^ = {D + rE)U, , 
where 14^ is a Lorentz like factor 



1 



1/2 



(44) 
(45) 

(46) 

(47) 

(48) 



and r is an adiabatic index for the equation of state, [Eq. 

In terms of these state variables, the hydrodynamic 
equations are as follows: The equation for the conserva- 
tion of baryon number takes the form. 



dD d log ( 



(49) 



The equation for internal energy conservation becomes, 

b^EV^) 



dE d log ( 



-P 



1 d 
d 



dt (j)^ dxi ' 



(50) 



Momentum conservation takes the form, 



dt 



-6S, 



91og( 
dt 



1 d 



dP 



-r.i^irs^v=)-^i^ (51) 



5, 



dx^ 
dx'- 



-WiD + TE)^^-aWiD + TE)§^^ , 

where x is the radiation reaction potential, which is 
described in § IIE . Note that the repeated occurance of 
the 4)^ factors simply account for proper volume factors 
(proper volume = (j>^{dx^)^). This is the reason for the 
choice of conformal scalings introduced in equations ( p^ ) 
and (p^. That is, we preserve D E, and Si when (f) is 
changed. 
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Our routines for evolving the hydrodynamics have been 
previously very well tested at the special relativistic level 
in ||3^,^,^,|9[^ where the hydrodynamics method de- 
scribed here was used to simulate relativistic heavy ion 
calculations. In that work, shock- wave solutions were 
compared for both decelerating and accelerating shocks. 
For the former case, the numerical results were accurate 
to better than 1% over the range of special relativistic 
7-factors from 1 to 10, i.e. Ethermai =0 — 10 pc^. For 
the case of shocks accelerating matter, errors increased 
to 1% for 7 > 2. 

A shock tube calculation was made using a constant 
adiabatic index of F = 2. The density ratio was 100:1, 
and the initial thermal energy density of the dense matter 
was equal to the initial baryonic density. The compres- 
sion ratio of the shocked bow density material was 8% 
too high. The rest of the density profile was accurate to 
better than 1%. Relativistic shock tube solutions were 
also calculated to test the accuracy of the rarefactions. 
Again in the range of interest {Ethermai ~ pc^) the over- 
all agreement of numerical results with exact solutions 
was of order 1%. 

On the basis of these results, we anticipate all shocks 
occurring in the neutron star coalescence will be treated 
with sufficient accuracy. Numerical errors are not 
Lorentz invariant, but tests for invariance have shown 
that our numerical methods are reliable for changes of 
frame by a factor of a few in rapidity. In the present 
work we have extended the hydrodynamics to curved 
space. However this extension is straightforward and we 
do not anticipate any new instabilities or unacceptable 
inaccuracies to be introduced thereby. For the calcula- 
tions described here there is little fluid motion and no 
strong shocks in our rotating frame. 



D. Equation of State 

For the orbital calculations presented here we use the 
zero temperature, zero neutrino chemical potential equa- 
tion of state from the supernova numerical model of 
Mayle and Wilson |^,^ . While the orbital calculations 
of concern in this paper should only involve zero temper- 
ature, there is some small shock heating of the stars as 
they adjust to changing conditions on the grid. Thus, 
we augment this equation of state with a thermal com- 
ponent (taken to behave as a F = 5/3 gas) in order to 
follow the dynamic evolution equations. Thus, we write 



P = Poip) + 7;Pie - eo{p)) 



(52) 



where Pq and eo are the zero-temperature pressures and 
energies. 

Table 1 gives the the zero temperature values of Pq, 
eo, and F vs. p. In figure 1 we present the bary- 
onic mass Mb, gravitational mass Mq, the stellar ra- 
dius in Schwarzschild (rg) and isotropic (r/) coordinates. 



the lapse function a, and the total mass energy density 
p{l -I- e) as a function of the central density of an iso- 
lated neutron star. From this figure it can be seen that 
this equation of state gives an upper limit to the gravi- 
tational mass of an isolated neutron star of Mq < 1-60 
Mq . [Note that this value supersedes the lower value pre- 
viously quoted in ]3^.] This limit roughly agrees with the 
upper limit of the smallest range of neutron star masses 
which overlaps all observational determinations. The fact 
that this upper limit is close to the typically observed 
neutron star mass Mq ~ 1.45 Mq has important conse- 
quences for the examples considered below. 

The three dimensional calculations reported on here 
have only about 15 zones in radius to represent each 
neutron star. As another test of the accuracy of our 3 
dimensional calculations, therefore, a hydrodynamic cal- 
culation was made of a single star using typical zoning in 
three dimensions. This calculation was compared with a 
one dimensional spherical hydrodynamic calculation with 
fine zoning. For the same baryonic mass, 1.60 Mq, the 
gravitational masses agreed to 2/3%, i.e. yielding a grav- 
itational mass of 1.45 Mq and 1.46 Mq for the 3D and 
ID calculations respectively. This we take as indicative 
of the accuracy of the calculated gravitational binding 
energy of the binary system as well. Figure 2 shows the 
density profile for a single isolated Mq = 1.45 Mq neu- 
tron star with the adopted equation of state.. 



E. Gravitational Radiation 

In general it is possible to express the emission of grav- 
itational radiation in terms of an "exact" expansion of 
multipole moments of the effective stress energy tensor, 
including corrections for the so-called "slow motion" ap- 
proximation . It is important to appreciate that these 
formulas apply to strong-field sources as well as to weak 
field sources |4^,Q| as long as the relevant components of 
the effective stress energy tensor can be identified. Since 
in the present paper, we are only concerned with orbital 
motion of equal mass binaries, the multipole expansions 
reduce to only a few nonzero terms. These we evaluate 
and test for convergence of the expansion. We summa- 
rize below the aspects of Q| which are relevant to our 
model. 

In any coordinate system (such as the one we are using 
here) in which the gravity waves far from the source can 
be characterized as linear metric perturbations propagat- 
ing on a flat background, the transverse traceless part of 
the metric perturbation characterizes the radiation com- 
pletely. This metric perturbation can be expressed [ p4[ in 
terms of the mass multipole (/'™) and current multipole 

Irn \ 



moments {S ) as 



EE 

1=2 m=—i 



-1 {I) rim 



/""(i-r)r 



E2,lm 
3 k 
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(53) 



where the superscript TT denotes the transverse traceless 
part of the metric perturbation and the notation (Oj'™ 
and ^''^ iS*'™ denotes the Ith time derivative of the respec- 
tive moments. 

From this, the general expression for energy loss is 



327r 



1=2 m=~l 



where the brackets denote averages over several wave- 
lengths. Angular momentum loss can similarly be writ- 
ten 

dt 32tt ^ ^ ^ ' 

1=2 m=-l 

+ ^(05';™*^(i+i)5'i™) ^ (55) 

where the expression Eq. (^5|) assumes an alignment of 
the angular momentum vector with the z axis. 

The radiation reaction potential x for Eq. ( ^7| ) can be 
written 

oo / 
1=2 m=-l 

Our problem then reduces to the identification of the 
relevant mass and current moments in our coordinates. 
For an asymptotically Minkowski coordinate system one 
can define a quantity 



(57) 



where g is the determinant of the metric and 7]°^^ is the 
Minkowski metric tensor. If h"^ satisfies the de Donder 
gauge condition 



(58) 



then the Einstein field equations take the form 

ah^f^ = -16TTT"f^ , (59) 

where t"^ is the " effective stress-energy tensor" [M . 

As long as the de Donder condition is valid, Equa- 
tion ( p9| ) can be inverted (using the fiat-space outgoing 
Green's function) and the Green's function expanded in 
terms of vacuum basis functions. The resultant expres- 
sion can then be reduced to provide expansions for 
the desired mass and current moments. 



IGtt 



(Z + !)(/ + 2) 



(2/ + l)!!V 2(1-1)1 

^ 167T 

k=0 ^ ' 



1/2 



\2k 



ipqi 



J+2k 



{2l + 2k+l)f (Z+l)(Z + 2) Y^^^2l^2,ln» 



2(fc + l) V2(2/-l)(2;-hl) 



pq 



f i{l-l){l + 2) V'^ 



2 l.lm^ 



\{2l - 1){21 + 3) J 

~ ^) \ ^ rp2 l+2,lm* 

2{2l + l){2l + 3)J 



2k 



21 + 2k + 3 
(fx , (60) 



and 
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where the Y * are the usual spherical harmonics, and 
Tpq'^"^* are the pure-orbital tensor harmonics as defined 
in 0. The first integral in Eqs. (|^) and (|^) is the 
usual spherical harmonic expansion. At the I = 2 level, 
Eq. ( |60| ) reduces to the well known quadrupole a ppr ox- 
imation. The second integral in Eqs. ( |60| ) and (^) is 
the correction to the slow motion approximation, which 
is non negligible in the present application, i.e. v/c^Q.l. 

To evaluate the time derivatives of the mass and cur- 
rent multipole moments we make use of the rotation 
properties of spherical tensors whereby, rotations can be 
generated in terms of the Wigner D matrices, 



jhn 7-)/ jln 



elm T^l nlm 



(62) 



where /q™ and 5*0™ are evaluated in the rotating frame 
For stable orbits (neglecting gravitational radiation) and 
hydrostatic stars, these are time independent quantities. 

The main contribution to the time derivatives is that 
due to orbital motion. Evaluation of the orbital motion 
reduces to derivatives of the Dl^^^, which for our coordi- 
nates have a simple ^ cos{mLut) dependance. 

The problem with evaluating Equations (^^ and ( |6l| ) 
is that the multipole moments are only defined in the 
de Donder gauge and not for our conformally flat coor- 
dinates. Furthermore, even if the transformation to our 
coordinates were straight forward (which it is not) the 
effective stress energy tensor would not be known. 

Fortunately, however, a transformation to de Donder 
coordinates is not necessary. It is only necessary that the 
moments of the metric coordinates be defined in a coordi- 
nate system which, like a de Donder coordinate system, is 
asymptotically Cartesian and mass centered (ACMC). In 
UM it is proven that in such coordinate systems and the 
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covariant metric components are time independent and 
expandable into a spherical harmonic (1/r) structure in 
terms of the same moments (i.e. Eqs. |6^ and|6|) relevant 
to the radiation field. Furthermore, these multipole mo- 
ments are invariant under transformations between two 
ACMC coordinate systems. From these expansions we 
can deduce the source for the slow-motion moments to 
be used in the equations for the radiation field, (|5^^5|) 
For example, the spatial three-metric must obey 




^ 1 



/=0 



{2i-iy.\f 2{i-i)i 



(/ + l)(/ + 2) 



1/2 



J2 jlmylrn + _ 1 p^le) + ■ ■ ■ + (0 polc) 



(63) 



On the other hand, our spatial three metric (Eq. ^ can 
also be expanded as the fourth power of a multipole ex- 
pansion of the flat-space Poisson equation for (p (Eq. f 9). 



where 



E 



E (21 

1=0 m=~l ^ 



An 



/■my;m^-(i+l) 



I \rlrn* 



(64) 



(65) 



and pi is the source term for (Eq. E0[) . If we collect the 
dominant linear terms in Equations ^6^) and (|6^ ) accord- 
ing to the recipe given in then we can identify the 
relation between the source pi for the conformal factor 
elliptic equation ( p^ ) and the mass multipole moments. 



327r 



(/ + 1)(Z + 2)V/' 



(2; + l)!!V 2(^-1)/ 



(66) 



This identification also reduces to the correct Newtonian 
limit. As can be seen from Eq. (^0|), pi p/2, where 
p is the Newtonian matter density, so that tqq — > p as 
required. 

The contribution from the current moments is ex- 
pected to be small as is the slow-motion correction. 
Therefore, we are mainly concerned with estimating the 
magnitude of those contributions. To the accuracy de- 
sired, we identify the source for the current moments S'''" 
and the slow-motion corrections with the Newtonian-like 
counterparts, i.e. we set Toj = Tqj, tij = Tij. We com- 
pute terms out to uj^'^, which includes mass multipoles 
out to Z = 4, current multipoles out to I = 3 and the 
leading correction for the slow motion correction. 

In Table |l| we summarize the relative contributions of 
various moments to the energy and angular momentum 
loss rates. As expected, the quadrupole term dominates. 
The next largest term is the slow motion correction which 
contributes only a few percent to the gravitational radi- 
ation and tends to decrease the loss rate. 



F. Numerical Methods 

The elliptic equations for the field and differential evo- 
lution equations for the hydrodynamic variables were 
finite-differenced in a Cartesian grid. The intrinsic state 
variables, D,W,E,T,a,(j>, are treated as zone centered 
quantities, while the four velocity Ui and momentum 
densities Si are node centered. The shift vector /3* and 
the three-velocity V are face centered. After finite dif- 
ferencing, the elliptic equations are reduced to a matrix 
equation. 



M ■ X 



(67) 



where M is a sparse matrix, x is a vector representing the 
relevant field variable at each zone, and b is derived from 
the source terms. This equation can then be solved using 
any one of a number of fast matrix inversion techniques. 

When we solve the elliptic equation for 0, the coordi- 
nate density D is adjusted so as to preserve the confor- 
mal scalings, Eqs. (0) and (|l|). That is, D = (jj^D is 
kept constant, which preserves baryon number. Also, the 
coordinate energy density is changed to preserve cj)^^ E 
and the momentum density is changed to preserve 4>^ , 
which maintains the entropy. 



1. Extracting Physical Observables 

The gravitational mass we obtain from the asymptotic 
behavior of ^ 1 + {GM/2r) (cf Eq. |6|). The angular 
momentum is more difficult to define. We estimate this 
from an integral over the space time components of the 
stress energy tensor 41 1 neglecting angular momentum 
in the radiation field. 



T'°x^ ~P°x' ]dV 



(68) 



Aligning the z axis with the angular momentum vector 
then gives, 

J = / {xSy ~ yS-')dV . (69) 



2. Boundary Conditions 

As noted above, our choices for the metric and slicing 
condition lead to a form for the Hamiltonian and momen- 
tum constraints in terms of flat-space elliptic equations, 
i.e. Eqs. (|l9|), (27), and (^), for the metric variables 
4>, (acj)) and A solution to these elliptic equations, 
however, requires that we specify values for (f), (acf)) and 
along the outer boundaries of the grid. For a Pois- 
son like equation, the field variables could be specified 
by integrating the source function over the interior, e.g. 
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(70) 



However, the evaluation of this integral for each point 
along the boundaries is computationally slow. In prin- 
ciple, an expansion of the source function in spherical 
harmonics Y^"^{9, (f>) could be applied to obtain the field 
variables along the boundaries, e.g. 



= 1 



EE 



An 
21 + 1 



1) qlrriYl™- ^ 



,0) . (71) 



However, the convergence of a spherical harmonic ex- 
pansion for a source dominated by two separate nearly 
spherical distributions is quite slow. In order to accom- 
modate these two features, we employ a combination of 
them which is numerically efficient even at distances rel- 
atively close to the neutron stars. 

In the equations for </> and (acj)) the boundary values 
are dominated by contributions from the effective point 
source potential from each star. Our method of speci- 
fying the boundary for (j) and (q;0), therefore, is to first 
make a best fit to the source density, pi or p2, of each 
star with a truncated spherical gaussian profile located 
at the source-density center of mass in each half of the 
grid. The boundary values on the computational grid 
then begin with the sum of the two point-mass contribu- 
tions from the truncated spherical gaussian profiles for 
each star. This provides a simple analytic contribution 
around the boundary for the bulk of the source. 

These gaussian density profiles are then subtracted 
from the source density to yield a residual density. An 
expansion in spherical harmonics up to Z = 4 (Eq. |7l| ) 
is then utilized to compute the contribution from the 
residual source function over the grid. For stars well sep- 
arated on the grid, this residual source typically accounts 
for only a few percent of the total boundary value. Also, 
the spherical harmonic expansion for the residuals con- 
verges faster than for the unsubtracted source function. 
Hence, a truncation of the expansion to Z = 4 is suffi- 
ciently accurate. 

A gaussian source profile turns out to be an excel- 
lent approximation in the early stages before the neutron 
stars begin to coalesce. In future calculations in which 
the stars will be followed until they merge, however, they 
will more closely represent a single source function. At 
some point in the calculation it will become expedient, 
therefore, to apply the spherical harmonic expansion di- 
rectly to the unsubtracted source function. 

We note that the expansion of the three metric (Eq. 
|63|) requires that the asymptotic form for obey. 



1 



2r 



Similarly, from the ACMC expansion for t/oo, 
lapse function must approach. 



1 - 



ma 



(72) 
] the 

(73) 



in order that our time coordinate become proper time as 
r — > cxD. The Poisson equation ( ^ ) for (g^) can also be 
expanded in spherical harmonics (e.g. Eq. ^l|) yielding 



{a4>) {aq 



(74) 



where ttiq^ is the volume integral over twice the source 
P2- Since in general m^^ ma, we choose the boundary 
condition, 



{a(f))c 



ma 



(75) 



to guarantee that Eq. ([73|) is satisfied. For the systems 
studied here, (a(/))oo ~ 0.98. 

In our computation of the boundary conditions, we 
impose a spherical cut off in the matter distributions 
at a radius equal to the largest sphere that fits within 
our cubical grid. This avoids the possibility of a spuri- 
ous hexadecapole moment associated with the cubic grid 
employed in the calculation. For matter terms this is a 
reasonable truncation for the calculations presented here, 
since only a negligible amount of matter appears beyond 
the surface of the neutron stars. However, the KijK'^^ 



terms in Eqs. (|l9|) and (27) contribute beyond the mat- 
ter boun dary . Also, the shift vector elliptic equations, 
( ^ ) and ( |3q ) , involve a source which extends beyond the 
source boundary. 

Regarding the KijK^^ terms we note that these terms 
are small. For example, the contribution to the gravita- 
tional mass from an integration over the interior source 
function is only ~ 0.0001 Mq. Furthermore, the asymp- 
totic form for KijK^^ should decay as l/r^. Assuming 
this form, we estimate that the exterior contribution from 
the KijK'^^ term is ^10~^ Mq and can therefore be ne- 
glected in the examples considered in this paper. 

Regarding the solution for the shift vector (Eqs. ^and 
|3^ ) we note that V • /3 is small and changes sign across the 
grid. This means that the variable x asymptotically goes 
to zero. Hence, we impose x = along the boundary for 
Eq. |3^. A solution for requires that we specify the 
boundary condition for the "drag" component G*. For 
this we note that G" behaves as an angular momentum 
density and should scale along the boundary as 



AyJ 



AxJ 



(76) 



III. ORBIT CALCULATIONS 

It is a nontrivial endeavor to find initial configura- 
tions for the two neutron stars prior to coalescence. Our 
method consists of placing two neutron stars on the grid 
with a rotational velocity sufficient to keep them in orbit 
and an initial "guess" density profile from a solution to 
the Tolman-Oppenheimer-Volkoff like equation for two 
single neutron stars in our isotropic coordinates. The 
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conversion from single star solution to a binary solution is 
achieved by allowing the stars to relax to an equilibrium 
configuration on the grid. That is, the field equations are 
then solved and the hydrodynamics evolved (without the 
radiation reaction potential and with viscous damping of 
the fluid motion) until equilibrium is achieved. For the 
examples to be presented below, we follow the time evo- 
lution of the system with constant angular momentum 
until it has settled down. As the stars settle down the 
damping is slowly removed. 



IV. RESULTS 

In this paper we are presenting principally the method 
of solving the field equations and hydrodynamics for bi- 
nary neutron stars. As examples, we also discuss below 
three illustrative calculations made at selected values of 
the orbital angular momentum with no radiation damp- 
ing of the orbits. Hi ghli ghts of these results have been 
presented previously ||35| . Here, we supply more details 
of the application of the method. 

In these calculations the neutron stars are taken to be 
of equal mass. The baryonic mass was selected so that in 
isolation each star has a gravitational mass of 1.45 Mq. 
Although the calculations presented here ignore radia- 
tion damping, during most of the evolution the radiation 
damping is small. Therefore, the stars should follow a se- 
quence of quasi-equilibrium configurations which closely 
match the equilibria computed here. These equilibria can 
be analyzed to obtain the rate of energy and momentum 
loss. Ultimately, the implied orbit decay could be used 
to infer the approximate time evolution through this se- 
quence of quasi-equilibrium orbits. 

We have placed the stars at various separation dis- 
tances on the grid and only run the calculation long 
enough to check whether the system obtains what would 
be a stable orbit in the absence of gravitational radiation. 
In total for the three orbit calculations presented below 
we have followed the stars through more than 20 revolu- 
tions (with roughly two hours of CRAY/YMP CPU time 
per orbit) to insure that the orbits have had time to set- 
tle down. We have utilized a grid of 100 x 25 x 25 zones 
for the matter and 100 x 50 x 50 for the field variables. 
We make use of reflection symmetry in the orbital plane. 
Also, since here we study equal-mass binaries, we exploit 
reflection inversion symmetry through the axis joining 
the centers of the two stars. In effect, then this calcula- 
tion is equivalent to a three-space grid of 10^ zones. 

Initial conditions for two 1.60 Mq baryonic mass neu- 
tron stars were obtained from the isotropic-coordinate 
form of the Tolman-Oppenheimer-Volkoff hydrostatic 
equilibrium. Two single star solutions were then placed 
on the grid and the hydrodynamic and fleld equa- 
tions evolved with viscous damping until a new quasi- 
equilibrium configuration was obtained. In practice this 
need only be done once for a selected initial angular mo- 



mentum. Subsequent orbits are then calculated by per- 
turbing the angular momentum and allowing the system 
to relax to a new stable configuration (if one exists). 

The first calculation was made with an orbital angular 
momentum of 2.2 x 10^^ cm^. The stars settled down 
into what appeared at first as a stable orbit, but later 
(less than one complete orbit) the stars began to slowly 
spiral in. For this system the angular momentum was not 
enough to support the orbit. Parameters characterizing 
this binary at the final time slice calculated are given in 
[Note that this table supersedes the table in 



Table [I] 



[ |35[ where lower values were quoted for the central densi- 
ties.] The stars were followed to a coordinate separation 
di « 34 km which corresponds to a ratio of proper dis- 
tance to total isolated mass of the system m = 2Mq of 
dp/m = 9.2. Where Mq is the gravitational mass of a 
single isolated neutron star. By this time it could be con- 
cluded that no stable orbit would result. Note that the 
minimum coordinate light speed a/cf)'^ is only 0.23 in this 
case. 

Figure ^ shows x, y contours in the z = plane for 
the hydrodynamic density D and the metric variables 
a, and 0^. Countours are drawn for the final time slice 
calculated for this angular momentum and the other two 
cases studied. 

We note that even at the last time calculated the stars 
are still quite far apart i.e. the ratio of coordinate sep- 
aration distance dj to their single-star radius, dj /rj ^4. 
Even in a Newtonian hydrodynamic calculation W5| , one 
would expect at most a few percent tidal distortion at 
this distance. The distortion is made even less, however, 
by the strong relativistic gravitational field (depicted in 
a and 0^ in Fig. ^) around the stars. We also note, that 
the central density increased continuously as the stars 
spiraled in. The orbit decay time, however is shorter 
than the collapse time. Thus, it seems likely that nei- 
ther the stars nor the orbit are stable for this system as 



summarized at the bottom of Table III 



The next calculation was made with an angular mo- 
mentum of 2.3 X 10^^ cm^. The orbit now appeared sta- 



ble as summarized at the bottom of Table III. However 



after about 1 to 2 revolutions the central densities were 
noticed to be rising. By the end of the calculation the 
central baryonic densities had continuously risen to about 
2.7 X 10^^ g cm~^ (« 10 times nuclear matter density) 
which is near the maximum density for a stable neutron 
star (cf. Fig. |l|). It appears that neutron stars of this 
mass range and the adopted equation of state may con- 
tinue to collapse as long as the released gravitational en- 
ergy can be dissipated. For this orbit the stars are at 
a separation distance of dp/m — 9.5, far from merging. 
However, the nonlinearities in the gravitational field have 
pushed the stars over the critical density for forming a 
black hole. By the time the calculation was ended, the 
minimum a had diminished to 0.379 and 0^ risen to 2.05 
corresponding to a minimum light speed of 0.18. Fig- 
ure ^ shows the stars as somewhat more compact objects 
which are continuing to collapse. 
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The final calculation was made with the angular mo- 
mentum increased to 2.7 x 10^^ cm^. As can be seen 
in Table III and figure ||, the stars at this separation 
dp/m = 12.4 appear both stable and in a stable orbit. 
Note, that even for this distance, the gravitational field 
(a and 0) remains strong. 

Figure ^ shows vector fields for Ux,Uy, Vx,Vy, and 
(ix,f3y for the calculation with J — 2.3 x 10^^ cm^. Fig- 
ures for the other two runs look quite similar. The max- 



imum value of the covariant velocity is \Ux,Uy 



0.9. 



The contravariant three velocities show the net fluid mo- 
tion after subtraction of the orbital motion. The max- 
imum magnitude is |T4,V^| = 0.07. We note that even 
in the corotating frame, Vx,Vy exhibits some rotation. 
Although the stars are initially corotating, as the stars 
relax to their equilibrium orbit, they acquire net rota- 
tion relative to the orbital motion. The net fluid motion 
is still much greater than the collapse velocity which is 
more than an order of magnitude slower than the maxi- 
mum vector shown here. A peculiar aspect of the velocity 
fleld in the rotating frame is that the fluid appears to be 
circulating about a vorticity which does not coincide with 
the central extrema in density or fleld variables. It also 
appears that the parts of the stars most distant from the 
companion rotate in an opposite sense, generating a sec- 
ond vorticity. The shift vectors are also plotted with the 
rotation subtracted. What remains is frame drag around 
the star due to its rotation which has a maximum ampli- 
tude of 0.015. It is also slightly offset from the centers of 
the stars. 



A. Analysis of the Collapse Instability 

Having observed the collapse instability numerically in 
the neutron star binary, one of course would like to have 
at least a qualitative understanding of the source of this 
deviation from Newtonian intuition. Here we present a 
heuristic explanation of the observed increase in density 
as the stars approach each other. We trace this increase 
to the effects of the Lorentz-like factor — 1. This fac- 
tor accounts for the speciflc kinetic energy of the orbital 
motion of the stars. Its effect is to increase the effective 
source strength. 

From Eq. (|l^ the Hamiltonian density pH has a term 
(\Y'2 _ l)(p -f- peT) which enters into the source term pi 
for (f). Similarly,the source for the Poisson equation for 
(acj)) (cf. Eq. |^ has a term {W^ - 1) x 3(p -h peT). 
Thus, the source terms for both (f) and a will increase as 
the separation distance decreases and W exceeds unity. 
A stronger source term will imply larger values for both 
(j) and a at the centers of the stars and therefore steeper 
gradients of these quantities as one moves outward from 
the center of each star. 

In isotropic coordinates, the general relativistic condi- 
tion of hydrostatic equilibrium for each star can be in- 
ferred from the dominant terms in the momentum equa- 



tion (51), 
dP 



(9 log a 



9 log a 
dx' 
d log 4 



dx^ 



dx^ 



(VF' - 1) 



(77) 



where we have ignored the centrifugal term, Sj{d(3^/dx^). 
From this we see that the effective gravitational force 
(right-hand side of Eq. 77) increases both because 
{W^ — 1) exceeds unity and because the gradients of a 
and </) are more steep as increases. A further in- 
crease of binding arises from the K^^Kij terms in the 
fleld sources, but these terms are much smaller than the 
— 1 contributions. 
In our rotating coordinate system, the fluid three ve- 
locities are nearly zero. Hence, from Eq. (E^) we have 



1 



(78) 



where R is the distance from the center of mass. Along 
the line between centers in the J = 2.3 x 10^^ cm^ model, 
an effective velocity of [luR^P' /a) = 0.28 is obtained. In 
the stars the average value is {W'^ — 1) ^ 5-10%. Includ- 
ing both the a and 4> terms in Eq. (77), we then estimate 
that the effective hydrostatic gravitational force on the 
stars is increased by 10-20% over that of stationary non- 
orbiting stars for which {W^ — 1) = 0. This increased 
gravitational attraction is sufScient to induce the calcu- 
lated increase of central densities as the stars approach. 
An important ingredient in this analysis, however, is that 
we have assumed zero thermal energy as the stars col- 
lapse, that is, we have assumed that neutrino emission 
is efficient enough to radiate away the released gravita- 
tional binding energy so that the stars remain effectively 
cold. 



B. Comparison with Post-Newtonian Results 

Much of the early evolution of a neutron star bi- 
nary should be describable with post Newtonian tech- 
niques. However, one desires an understanding of where 
in the evolution the post-Newtonian approximation di- 
verges from a fully relativistic treatment. For this reason 
it is of interest to compare the results here with those 
obtained by a post-Newtonian treatment. We caution, 
however, that such a comparison is ambiguous. The two 
formalisms invoke different gauge choices. Hence, param- 
eters can have different meanings. 

Our intermediate orbit (J = 2.3 x 10^^ cm^) appears 
to be on the verge of the transition from steady inspiral 
to unstable plunge. Therefore, it is convenient to com- 
pare these results with the (post)^/^-Newtonian analysis 
of this transition in an equal mass binary as given in ref. 
p3| . In that paper a search was made for the inner most 
stable circular orbit in the absence of radiation reaction 
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terms in the equations of motion. This is analogous to 
the calculations performed here which also have analyzed 
orbit stability in the absence of radiation reaction. 

In the (post)^/^-Newtonian equations of motion, a cir- 
cular orbit is derived by setting time derivatives of the 
separation, angular frequency, and the radial acceleration 
to zero. This leads to the circular orbit condition lESl, 



ujI = mAo/dl 



(79) 



where ujq is the circular orbit frequency and m = 2Mq, 
dh is the separation in harmonic coordinates, and Aq is 
a relative acceleration parameter which for equal mass 
stars becomes, 



A 3 771 

An = 1 

2dh 
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8 dh 
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(80) 



Equations (|7|) and (|o|) can be solved to find the orbit 
angular frequency as a function of harmonic separation 
dh- The gravity wave frequency is then twice the orbit 
frequency, / = wq/tt. 

Figure ^ shows the circular orbit post-Newtonian grav- 
ity wave frequency vs. separation distance compared 
with the present numerical calculations. A striking fea- 
ture of the present work is that as the stars approach 
one another, the frequency becomes nearly independent 
of the separation distance until the orbit becomes unsta- 
ble to plunge at a relatively large separation. 

The main parameter characterizing the last stable or- 
bit in the post-Newtonian calculation is the ratio of coor- 
dinate separation to total mass (in isolation) dh/m. The 
analogous quantity in our nonperturbative calculation is 
proper separation to gravitational mass, dp/m. The sep- 
aration corresponding to the last stable orbit in the post- 
Newtonian analysis does not occur until the stars have 
approached 6.03 to. For Mq = 1.45M0 stars this would 
correspond to a separation distance of about 26 km. In 
the results reported here, however, the last stable orbit 
occurs somewhere just below 9.4 itiq at a proper separa- 
tion distance of dp !=s 40 km. 

The fact that we observe the last stable orbit to oc- 
cur at larger separation is consistent with the results of 
Wex and Schafer , who found that higher-order post- 
Newtonian terms beyond those considered in were 
significant. They found an unstable orbit at 7.50 m with 
large possible uncertainty. A larger separation is also 
consistent with the numerical initial-data calculation for 
two black holes by Cook |^^ . In that paper a minimum 
proper separation between horizon's of 4.88 to was found. 
If the black hole result were naively applied to the curved 
space around the two neutron stars by adding this dis- 
tance to their Schwarzschild radii, this would correspond 
to a separation of 8.88 to between centers. 

The differences in J/2MqIi^ and angular velocity 
also refiect this larger distance. Note, however, that even 
though the stars are much farther apart, the relativistic 



treatment gives a stronger gravitational binding energy 
for this system. However, our binding energy includes 
the increased self binding of the individual stars. 

Given that the separation distance is greater and that 
Lu ~ r^/^, we would observe a frequency which is a fac- 
tor of (40/26)"^/^ « 2 slower than the post Newtonian 
estimates if separation distance were the only relevant 
effect. However, our angular frequency is about a fac- 
tor of two less than the post-Newtonian value even at 
the same separation distance. We can understand this 
nonlinear effect of general relativity heuristically. This 
difference in orbit frequency should trace to the balance 
between the gravitational and centripetal forces. From 
the momentum equation (51) with W = 1 we can write 
the dominant terms as; 



(D + TE) 



9 log a dp- 
:■ — = — r 



(81) 



The left-hand side is the analog of the Newtonian grav- 
itational force. The right-hand side is the analog of the 
centripetal force. 

Along the y = axis we have (3^ « lux and Sy = 
(D + TE)Uy. Thus, we may write 



9 log a 



(82) 



Taking the fluid three velocities in the rotating frame to 
be zero, Vi = 0, then from equation (^6|), 



uxcf) 



yjl — Uj'^(f>'^X^ 



Hence, 



9 log a 
dx 



■ IjP'X- 



(83) 



(84) 



From this we see that the analog of the Newto- 
nian centripetal term is enhanced by a factor of 
(\)^ I (p? -y/l — w^^^'x^) . The average of this quantity along 
a line joining the centers of the stars is « 7.0. Thus, a 
much smaller value for the orbital frequency w provides 
sufficient centripetal force to maintain a stable orbit. 



V. CONCLUSION 

We have summarized a method to study the relativis- 
tic evolution of a binary neutron star system. We have 
illustrated the method by following the evolution of a 
close binary through several orbits for different angular 
momenta in the absence of radiation reaction. These re- 
sults show two new results which to our knowledge have 
not been reported previously. 

One significant aspect of these calculations is that the 
binary orbit becomes unstable at a much larger separa- 
tion distance (a factor of « 1.6) than that derived from 
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(post)^/^-Newtonian analysis. This implies that the non- 
linear effects of gravity become important much earlier 
on, and that searches for gravity waves may observe the 
final merger to occur at a lower frequency than expected. 
This is important, since it places the coalescence closer 
to the maximum sensitivity range of the LIGO detector 
and others. However, we estimate little tidal distortion 
or hydrodynamic amplification of the the gravity wave 
signal. Our estimate of the gravity wave amplitude near 
the final orbit is /i « 3 x lO'^^ (at 100 Mpc), (see Table 
||for a summary of the gravitational wave information) . 

A second significant result of the present calculation 
is that the nonlinear effects of the fully relativistic grav- 
ity imply deep gravitational wells as two neutron stars 
approach each other. Indeed, the fields can become so 
strong that the stars are individually unstable to collapse 
into two black holes. Exactly when or if this instabil- 
ity occurs is of course dependent upon the equation of 
state. However, for the realistic equation of state em- 
ployed here, this collapse is observed to occur while the 
stars are still in a quasi stable orbit. This suggests that 
there could be many orbits after black-hole formation 
until the stars actually merge as two black holes. If cor- 
rect, this result will have a significant impact on future 
studies of binary neutron star mergers and renders the 
two-black-hole coalescence much more important. Fur- 
thermore, the possibility of a collapse many orbit periods 
before coalescence may have observational consequences 
not only for gravity wave detectors, but in electromag- 
netic (optical, radio, x-ray, and 7-ray) signals as well. 
We note that even for our unstable inner orbit, the spe- 
cific angular momentum a = J/Mq w 1.3 Is greater than 
unity, implying that more angular momentum must be 
lost before merger can occur, although a transient black 
hole may be able to form with a > 1. 

The importance of the possibility of premerger stellar 
collapse is dependent on the equation of state and the 
distribution of neutron star masses. As stated earlier, 
the EOS we use is the same one used by Mayle & Wil- 
son [ ^Tp^ ] in their supernova model. Calculations made 
with this EOS for a model of supernova 1987A give an 
explosion energy of 1.5 x 10^^ ergs, consistent with ob- 
servation. Also, the neutrino spectra and time of neu- 
trino emission are in good agreement with the 1MB and 
Kamiokande neutrino detections . These models also 
give a good reproduction of heavy element nucleosynthe- 
sis in the baryon wind from the proto-neutron star . 
An important point is that with a stiffer EOS (which 
would allow a higher mass neutron star), Mayle & Wil- 
son were not able to obtain satisfactory results. 

Regarding, the up per mass limit to neutron stars from 
observations, Finn ||55| has recently analyzed the ob- 
served masses of neutron stars and has assigned a lower 
limit of 1.15 to 1.35 Mq and an upper limit of 1.44 to 
1.50 Mq at the la (68%) confidence level. At the 2a 
95%) confidence level the upper limit only increases to 
1.43 to 1.64 Mq. In an independent approach, Bethe 
and Brown [p6| have recently argued from nucleosynthe- 



sis constraints that the maximum neutron star mass is 
1.56 M0. They also point out that if kaon condensation 
is taken into account the critical mass may only be 1.50 
Mq. As seen in Fig. |l|, our upper mass limit is 1.60 
M0 and the mass of the stars in our sample calculations 
is 1.45 Mq, quite consistent with the observational lim- 
its. If the maximum observed stellar mass were as low 
as the la upper limit, i.e. 1.50 Mq, it could be that al- 
most all neutron star binaries would precollapse before 
coalescence. 

The sample calculations presented here were made 
with a relatively coarse spatial zoning (see §11 D). We 
therefore, consider the present results to be only qualita- 
tively correct. They must be supplemented with more de- 
tailed numerical modeling. We are presently developing 
methods of improved numerical efficiency so that results 
will be of sufficient accuracy to allow for a quantitative 
interpretation of future observed neutron star binary col- 
lapse and coalescence. 

Note that for the above examples J/ujJ ss 10"'' justify- 
ing our neglect of gravitational radiation. Also, we point 
out that since this paper was submitted it has come to our 
attention that a study has been made of the validity 
of the conformally flat condition when applied to rapidly 
rotating isolated neutron stars. This is the simplest case 
for which the conformal approximation is different from 
the exact equations. This work is encouraging in that it 
has shown that the method works remarkably well. 

From the above discussion it is clear that further stud- 
ies are warranted, particularly an effort to better deter- 
mine the last stable orbit and the approach to this or- 
bit, as well as a systematic study of the sensitivity of 
the collapse instability to the neutron star EOS. Work 
along this line is currently in progress. There is also a 
need to investigate orbits at larger radii so that a reliable 
connection to calculations in the post-Newtonian regime 
can be made. Once this is done, one can combine the 
post-Newtonian and numerical (3-1-1) results to produce 
a template of expected gravity wave signals. These can 
then be used to better extract the gravity-wave signal 
from the noise. 
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FIG. 1. Various parameters characterizing isolated neutron 
stars with the adopted equation of state as a function of the 
central baryon density pc. Mb and Ma are the baryonic and 
gravitational masses, respectively, in units of Mq . The radius 
is given in both Schwarzschild coordinates rs and in isotropic 
coordinates ri in units of 10 km. Also shown are central 
values for the minimum lapse function a, and the total mass 
energy density ptot = p(l + e) (in units of 10^^ g cm~^). 



FIG. 2. Density profile as a function of radius (in isotropic 
coordinates ri) for an isolated neutron star with a gravita- 
tional mass of 1.45 M©. 



FIG. 3. Contours in the orbit plane for various quantities 
at the final time slice calculated for J = 2.2 x 10^^ cm^, 
J = 2.3 X 10" cm^, and J = 2.7 x lO" cm^ as labeled. 
The various figures are for; a) coordinate density D = pW; 
b) lapse function a; c) conformal factor cj)^. Note that the 
J = 2.3 X 10^^ cm'^ density contours are more compact illus- 
tratin g th e degree to which these stars have collapsed. See 



TABLE 1. Zero temperature equation of state 



agtr 

x 1.09 



Table y| for pmax, (jilaax, and amin- Since Wmax = 1.09, 



FIG. 4. Vectors in the x, y orbit plane around the star cen- 
tered at a; = —17 A km. This is at the final time calculated 

The various figures 



for the binary with J = 2.3 x lO" cm^ 



are for: a) four velocity , Uy (showing the overall orbit mo- 
tion); b) the rotation subtracted contravariant three- velocity 
, (showing the corotation and collapse of the fluid); and 
c) the rotation subtracted shift vectors Px,l3y (showing the 
relativistic frame drag around the star). 



FIG. 5. Comparison of gravity wave frequency / vs. proper 
separation distance dp from the present work (points) and 
post Newtonian circular orbit condition (line). The arrows 
show approximate location of the last stable circular orbit in 
the two schemes. 
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TABLE II. Contributions to energy and momentum loss 
from the orbit calculation with J = 2.7 x 10^^ ' — 
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TABLE III. Parameters characterizing the orbit calcula- 
tions at the final edit. Mg is just the total mass of the binary 
divided by 2. 
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